Structure of binary Bose-Einstein condensates 



Marek Trippenbach^'^, Krzysztof Goral^, Kazimierz Rzazewski'^, Boris Malomed'* and Y. B. Band^ 

^ Departments of Chemistry and Physics, Ben-Gurion University of the Negev, Beer-Sheva, Israel 84105 
^ Institute of Experimental Physics, Optics Division, Warsaw University, ul. Hoza 69, Warsaw 00-681, Poland 
Center for Theoretical Physics and College of Science, 
Polish Academy of Sciences, Al. Lotnikow 32/64, Warsaw 02-668, Poland 
^ Department of Interdisciplinary Studies, Faculty of Engineering, Tel-Aviv University, Tel-Aviv, Israel 69978 

Abstract 

We identify all possible classes of solutions for two-component Bose-Einstein condensates (BECs) 
within the Thomas-Fermi (TF) approximation, and check these results against numerical simulations 
of the coupled Gross-Pitaevskii equations (GPEs). We find that they can be divided into two general 
categories. The first class contains solutions with a region of overlap between the components. The 
. other class consists of non-overlapping wavefunctions, and contains also solutions that do not possess the 

' symmetry of the trap. The chemical potential and average energy can be found for both classes within 

, the TF approximation by solving a set of coupled algebraic equations representing the normalization 

conditions for each component. A ground state minimizing the energy (within both classes of the states) 
is found for a given set of parameters characterizing the scattering length and confining potential. In the 
, TF approximation, the ground state always shares the symmetry of the trap. However, a full numerical 

solution of the coupled GPEs, incorporating the kinetic energy of the BEG atoms, can sometimes select 
a broken-symmetry state as the ground state of the system. We also investigate effects of finite-range 
interactions on the structure of the ground state. 
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! 1 Introduction 
^ . 

^SJ , Phase transitions and coexistence of different phases in multi-component systems are of great importance 

OO ' to many areas of physics, chemistry and biology. An ideal system to study these phenomena is a multi- 

I component dilute atomic gas Bose-Einstein condensate (BEG) mixture at zero temperature, due to the 

simplicity of its theoretical description. The mean-field approximation provides an excellent description of 
these systems. Other multi-component systems cannot be understood as well as these BEC mixtures, because 
' their density is generally much higher, and the delta-function pseudopotential, used to describe interactions 

between atoms in BEC, is not appropriate for them. Instead, a true microscopic interaction potential must 
be employed to adequately describe such systems, hence modeling them is much harder. 
' O I Multi-component BECs have been extensively studied over the last few years ||l|-||l0||. These studies have 

been motivated by experimental work performed by the JILA and MIT [p^Pgroups. Many interesting 
^ I effects have been experimentally determined and theoretically predicted, including topological properties of 

the ground and excited states, phase transitions and symmetry breaking effects produced by a phase 
difference between components ||] , stability properties , Josephson-type oscillations [|| , four- wave mixing 
' [Q, and trapping of boson- fermion and fermion-fermion systems Nevertheless, many features of BEC 

^ I mixtures remain to be explored by theorists and experimentalists. 

A large variety of different species can be used to produce mixtures of condensed bosons. Mixtures of 
two different elements, or of different isotopes of the same element, or simply different hyperfine states of 
the same atom can be considered. Simulating experimental results for BEC binary mixtures requires 

knowledge of the scattering lengths of the atoms involved. To the extent that the values of the scattering 
lengths are known with insufficient accuracy, a full classification of different states is necessary within the 
range of possible values. This is also necessary in the context of tuning the scattering length, as can be done 
by means changing the external magnetic field near Feshbach resonances Ref. Varying the scattering 

length, one can study phase transitions to states that break the symmetry of the trapping potential. Such 
states are known in the literature, and they were observed in the JILA experiment pl . 

The classification of two-component condensates can also be carried out in systems with interconversion 
of components ("chemical reactions" between them), as in the case of an atom-molecule condensate, where 
the conserved quantity is the number of atoms plus twice the number of molecules (the number of atoms 
and number of molecules are not separately conserved). A mathematical model of the latter system can be 
formulated in terms of two coupled Gross-Pitaevskii equations (GPEs), which contain, in addition to the 
familiar cubic self- and cross-interaction nonlinear terms, quadratic terms that account for the "chemistry" 
(i.e., the interconversion) |]l6| . In particular, an interesting prediction of the model is that a "soliton" state. 
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i.e., a stationary self-supported condensate cloud (similar to "light bullets" in nonlinear optics p7|), may 
exist without any trapping potential present p6| . 

Many aspects of binary condensate mixtures have been treated in the literature, using various (mostly 
numerical) methods in order to predict results for various experimental setups (see, e.g., Refs. |], § and 
Q). Nevertheless, a general classification of all the ground-state solutions is not yet available. This is 
understandable in view of many control parameters present in the models (three scattering lengths, particle 
numbers for both components and characteristics of the trap). A very general and elegant, but not explicit, 
algorithm for determining ground-state shapes has been proposed by Ho and Shenoy jl] . Here we start with 
the same goal in mind, but also with the intention to provide a maximally straightforward and analytic set 
of expressions for the ground-state wavefunctions and energies of atomic gas BEC mixtures. 

We use a computational method simulating the evolution of two-component GPEs in imaginary-time 
in order to study phase separation of components in BEC mixtures. Results produced by this method 
are compared with analytical predictions based upon the Thomas-Fermi (TF) approximation applied to 
two-component GPEs. We analyze the changes in the structure of separated phase BEC mixtures with the 
variation of the s-wave scattering lengths and atom numbers. The changes can be predicted and understood 
using a simple TF approximation, which includes equations obtained from normalization conditions for both 
components. The TF picture is compared with solutions obtained using numerically simulated CPE evolution 
in imaginary time, which includes the kinetic energy of atoms neglected in the TF approximation. We find 
that, for the simple case of a spherically symmetric harmonic trapping potential, many spherically-symmetric 
phase-separated geometries are possible, depending on the ratios of the self- and cross- s-wave scattering 
lengths for atomic collisions. In the TF approximation, symmetry-broken phase-separated geometries (i.e., 
those whose symmetry is lower than that of the trap) are always energetically higher in energy than those 
with unbroken symmetry. Nevertheless, numerical simulations in imaginary-time show that a lower-symmetry 
state may be the lowest energy eigenstate, and thus determine a ground state of the system. Our method 
may be generalized to include a finite-range interaction between atoms. In the last section of the paper we 
study, by means of direct numerical simulations, how such interactions affect the geometry and shape of the 
ground state. 



2 Mean-field description of two- component Bose-Einstein mix- 
tures 

In the present work, we concentrate on stationary states of BEC mixtures, (not their dynamics), therefore 
we start with the time-independent coupled GPEs, written in the standard notation: 

(-Mi-^ + Vi(r) + [/n|V;i(r,i)p + [/i2|V;2(r,i)|2^ Vi(r,t) = , (1) 

(-^^2 - 1^ + V2{r) + C/i2|Vi(r,0l' + f/22|^2(r,t)p) M^^t) = . (2) 

Here /zi_2 are chemical potentials of the two species, and Vi.2(r) are two isotropic parabolic trapping poten- 
tials, i.e., 

l/,(r) = (m,/2)w|r^ j = l,2, (3) 

LUj being the corresponding frequencies of harmonic oscillations of a trapped particle . Further, Uij = 
{^.ttTt? / rriij) Uij are atom-atom interaction strengths, proportional to the s-wave scattering lengths an, 022, 
and ai2 for the 1-1-1, 2-1-2, and 1-1-2 collisions, respectively, where 1 and 2 numerate the components, 
and ruij = mi if i — j and rriij = 11111712/ {mi + 7712) if i 7^ j. For simplicity, in the numerical calculations 
presented here we take mi = TO2 = m, and assume that the magnetic moments of atoms belonging to the 
different components are equal, so that the corresponding trapping potentials are equal too, uji = uj2 = uj, 
but this condition as well as the spherical symmetry condition may be readily relaxed by means of rescaling 
variables. Furthermore, we assume that the scattering lengths are real, i.e., we assume that collisions are 
not lossy. We also assume that all the scattering lengths for both different and alike atoms are positive; 
otherwise, the classification of the possible states becomes very cumbersome. 

Our calculations were carried out, simulating the evolution in the time- dependent GPEs in imaginary time 
, so that to let the solution relax to the ground state. The computations used the split operator method 
with the fast Fourier transform, similar to that used in Ref. [|l8|. The chemical potentials are obtained by 
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computing the net energy (the sum of kinetic, potential, self- and cross- nonlinear mean-field energies) for 
each component. We have chosen the wavefunctions Tpi{r,t) and ip2{'r,t) to be normalized to the number of 
particles in each component, so that J |'0i(r, t)\'^(Pr = Ni. Choosing the symmetry of an initial configuration 
in the imaginary-time simulations, a solution "01,2(1") which minimizes the total energy 

£ = l d3r[/ii|0i(r,i)p+M2|^2(r,i)p-^|0i(r,i)|4~^|V-2(r,t)|4-C/i2|^/^i(r,t)p|^2(r,t)P] , (4) 

within the class of functions possessing this symmetry can be found. The ground state of the two- 
component Hamiltonian is the one with the smallest value of £ for different symmetry classes. 



3 Thomas-Fermi Approximation 

The TF approximation can be used to describe a zero temperature condensate in the cases when the trapping- 
potential and mean-field nonlinear terms in GPEs are attractive and repulsive respectively, and the number 
of atoms is large so that the mean-field energies are large compared to the kinetic energy. For the two- 
component system with overlapping wavefunctions. 



/ii-l/i(r) \^{Un Ui2 \ ( |V'i(r,t)P 
M2-1^2(r) j I U21 U22 J \ \M^,t)\^ 



(5) 



These equations can be solved as a linear system of equations for |'0i(r,t)p and |V'2(r, OP terms of the 
chemical potentials fii and 112 to obtain: 



|V'i(r)|2 ^ L^:^::^^ — ^-i^^^ — r:- — v-,-,- ■ ^ 

U11U22 — U12U12 

|^^(j.)|2 ^ [M2^ii - Mi^i2] - [Uii - U12] {m/2)uj\^ 
U11U22 — U12U12 

The chemical potentials /ii and fj,2 are determined from the normalization conditions, Ni — J dP x |'0i(r, t)p, 
where D is the dimension. We shall plot examples for _D = 1, but our numerical method is valid for higher 
dimensions as well, hence we derive all the formulas for the general case. 

When phase separation occurs and there are regions in the physical space occupied by one component 
only, the TF approximation leads, instead of Eqs. (|5|), to the corresponding one-component TF equations. 
For example, if a phase with only species 1 exists in a particular region of space, an equation of the form 

|Vi(r,t)P = (/ii-yi(r))/C/n , (8) 

is to be used in this region. If another phase exists wherein the species 1 and 2 are mixed, Eqs. (^ are 
relevant for that region. The chemical potentials /ii and ^2 must be determined by setting the number of 
atoms of each type equal to the integral of the corresponding density over the whole space. 

Inspection of Eqs. (^ suggests then that all the solutions can be classified according the signs of three 
parameters: det J7 = U11U22 ~ (£^12)^, and aj = Ujj — Uu- In particular, q;i,2 determine signs of the 
curvature (coefficients in front of r^) of the effective quadratic potentials for the two components in Eqs. (^) 
and (^. Qualitatively different types of possible states with overlapping wave functions (i.e., disregarding 
regions where only one of the species is present) identified by the TF analysis are defined in Table 1. 

Let us focus on those cases when kinetic energy contribution does not change the general structure of 
the solution, but only generates narrow transient layers on the scale of the corresponding healing length, as 
in the single-component case when the TF approximation is valid (thus we consider large-size condensates. 
The analysis will include the TF configurations with and without the overlap of the two different condensate 
wavefunctions. As already mentioned, in the most cases both wavefunctions do not overlap everywhere, i.e., 
there is a region where only one wavefunction is different from zero. Consequently, search for the lowest- 
energy state of the mixture cannot rely solely on Eqs. (§|)-(0) obtained in the assumption that the overlapping 
takes place everywhere. 

Combining the cases represented in Table 1 and single-wavefunction solutions within the TF approxima- 
tion, we distinguish two general types of solutions: unseparated ones, having an overlap region where both 
wavefunctions are nonzero, and separated solutions which do not contain any overlap region. In the latter 
case, we shall see from the analysis of a full CP equation (with kinetic energy included) and also a nonlocal 
version of the two-species model, with a finite range of the interatomic interactions, that it is necessary to 
further distinguish between weak {U11U22 < ([^12)^) and strong {U11U22 ^ (^^12)^) separation 0. 



[MiC^ii - M2C/12] - [U22 - U12] {m/2)io^r^ 
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U11U22 — U12U12 


Uii - U12 


U22 — U12 




= detU 


= ai 


= a2 


Case 1 (Type A) 


positive 


positive 


positive 


Case 2 (Type B) 


negative 


negative 


negative 


Case 3 (Type A) 


positive 


negative 


positive 


Case 4 (Type A) 


positive 


positive 


negative 


Case 5 (Type B) 


negative 


positive 


negative 


Case 6 (Type B) 


negative 


negative 


positive 


Not possible 


positive 


negative 


negative 


Not possible 


negative 


positive 


positive 



Table 1: Classification of the Thomas- Fermi forms. 



3.1 Partially Overlapping Wavefunctions {U11U22 — (^12)^ > 0) 

In the case det C/ > (cases 1, 3 and 4 in Table 1), we have checked numerically that the minimum-energy 
solution is given by the wavefunctions of the form shown in Fig. |l|. Near the origin, both wavefunctions 
coexist up to the point where one of them vanishes. Past this point, one wavefunction vanishes, while the 
other one remains nonzero, following the single-component solution, Eq. (H) . Fig. |l| shows two different cases 
that are possible with the scenario described above. In Fig. |^a, the two effective trapping potentials have 
the same sign of their curvature, i.e., aia2 > 0, in the overlap region (this is case 1 in Table 1). Fig. [^b 
presents another situation, when the two effective potentials have opposite curvatures, aia2 < (these are 
cases 3 and 4 in Table 1). 

3.2 Separated Wavefunctions {UuU22 - (U^f < 0) 

This category is represented by cases 2, 5 and 6 from Table 1. The simplest configuration is that with one 
wavefunction being different from zero in the region around the origin and vanishing beyond a separation 
radius, R, while the second component surrounds the first one. In this case, we can express /ii and 112 as 
functions of R and minimize the net free energy £, in order to find the lowest eigenstate of this type. The 
normalization conditions for the wavefunctions of the two components give a set of relations between R, the 
chemical potentials fii and 112, and the number of atoms in each condensate: 

iVi = / d^r [fii - V{r)] /Un , (9) 
Ja 

i-Ro 

N2 = / d^r [^i2 - V{t)] /U22 ■ (10) 

JR. 

Here V{v) is the binding potential (^), and i?o is a outer radius at which the wavefunction of the second 
component vanishes in the FT approximation. We first consider the ID case and then show how these 
considerations can be generalized to two and three dimensions. 



3.2.1 One-Dimensional Case 



To find the value of the radius R minimizing £, we solve the set of the coupled equations (|^) and (10) for 
/ii and /i2. The first equation can be solved directly to yield /xi as a function of R. The second equation is 
more complicated - it can be solved analytically only in the ID and 2D cases, but not in 3D. In the ID (3D) 
case, one needs to solve a third- (fifth-) order algebraic equation to find IJL2 as a function of R. The system 
of equations (^ and (|l^) in ID reduces to 

2{^iiR-Knijj'^R^) ^UixNi (11) 

2 ( ^'^^^ - pi2R + l:muj^R^) ^ U22 N2 . (12) 

In the case under consideration here, without spatial overlap of components, the total energy is simply a 
sum of the average values of harmonic potential and half of the nonlinear term in the CP equation for each 
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component state: £ = ^ ^ T/^ilmw^x^ + C/M|V'iP/2|'0i >. If we substitute the direct expression for the 
wavefunction in the TF approximation we obtain: 



£ 



(13) 



and the chemical potentials can be found using Eqs. (]rT|)-(]T^). 
3.2.2 Generalization to Two and Three Dimensions 



A particularly simple form of these equations is obtained upon introducing a TF radius (i?TF)i 2 of the 
condensate in the corresponding dimension. The TF radius is defined as a radius of the single spherically- 
symmetric condensate obtained in the TF approximation. We can find an explicit dependence between the 
chemical potentials of the condensates and the separation radius of the two phases in 2D. In this case, we 
again define the TF radii, which in 2D are equal to: (-Rtfi)"' = SNiUu / {Trmuj'^) . In this case, we obtain the 
following set of equations for chemical potentials and total energy: 

,2 



Ml 



£ = 



4i?2 



2 



(RtfiT + R" 



(i?TF2) + 



1„6 
6 



1„6 
6 



p4 



p4 



(14) 
(15) 
(16) 



where r is defined as r = Rvmuj^ . 

In 3D we can define the TF radius is given by (i?TFi) 5 = {15Uu) / {2TTm(jj^) , and we can derive a set of 
equations which can be solved for the chemical potentials vs. the separation radius: 

Ml 



4(^)5/2 



10i?3 

3i?^- 



3i?'' + 2(i?TFl) 
. 2^2 



The energy is given in terms of the separation radius R by: 



i?3-2(i?TF2) =0. 



£ ^ 



1 



5(2;j2 
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(mti;2)^/2 



TFl 



-"-TF2 



(17) 
(18) 

(19) 



where r ~ RVtuuP. 



3.2.3 Symmetry Breaking Solutions 

Eigenstates of the binary-condensate system that break the symmetry of the trapping potential exist. In 
ID, a solution of this kind is given by TF parabolas that are stuck together. An example is shown in Fig. H. 



In this case one can derive equations for /xi and 112 in the same way as in Sec. 3.2.1, integrating the densities 
and substituting the result into the normalization conditions. The generalization to the dimensions higher 
than one may be only obtained if the separation surface (which reduces in ID to a single point) is simple. 
In ID, the corresponding coupled equations take the form 

(2u )^^'^ 1 

^ ^^'^ - fiiR + -muj^R^ = Uii Ni (20) 
3v mcij^ 6 

+ f^-2R Imc'R' = U22 N2 (21) 
3v mw^ D 

£ = (f^ - + l^imc.'fR^)/Un + (f^ + f^R ^{mu;rR')/U22 (22) 

We have checked numerically that this solution cannot give rise to a minimum of the free energy, hence, 
within the framework of the FT approximation, the ground state cannot be the one with broken symmetry. 
The difference in the energy between symmetric and asymmetric cases is usually very small making them 
almost degenerate. The degeneracy is exact in the limit when — > t/22 and is removed by the kinetic 
energy. 
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4 The role of kinetic energy 



In this section we present results of our studies of the contribution of the kinetic energy on the total energy 
and on the functional dependence of the ground state of the binary mixtures of the BEC. For a single 
condensate, in the regime of the validity of TF approximation, the kinetic energy creates a healing length, 
^, in the region where the condensate wavefunction tends to zero. This healing length is of the order of 
(87rna)^^/^ where a is a scattering length and n is an average density of the condensate. For the 

binary mixtures of BECs there are two length scales that can be defined in order to characterize two kinds 
of boundary regions. One is the ordinary healing length of single condensate and refers to the healing of 
the wavefunction outside region of the coupled wavefunctions. But for mixtures another region between the 
two components exist and a penetration depth, as a length scale over which two components overlap. 
The penetration depth is a function of det U. For det U < the lowest energy state consists of partially 
overlapping wavefunctions; hence the penetration depth is of order of the size of the condensate. With 
decreasing det U the penetration depth becomes smaller and goes to zero in the limit of strong repulsion as 
detU — > — oo. At the same time, the contribution of the kinetic energy to the total energy becomes more 
important, in spite of the shrinking overlap region. This situation is illustrated in Fig. ^. The energy of 
the lowest eigenstate is plotted in the symmetric and asymmetric classes as a function of det U / {U11U22) ■ 
Fig. 1^ Only one curve is plotted for det U > where the contribution of the kinetic energy is negligible and 
Thomas-Fermi approximation gives an excellent prediction for both the value of the ground state energy and 
its wavefunction. As the value of det U becomes negative, the lowest energy state within a TF approximation 
consists of two separated components and the energy does not depend on U12 and is plotted as a horizontal 
dashed line in Fig. |3[ det U < the contribution of the kinetic energy is substantial and it is larger for 
the symmetric case which has two interfaces between phases. The asymmetric configuration has only one 
interface in the lowest energy state. Hence, the ground state looses symmetry of the trap. 

In order to search for stable symmetry-breaking solutions we kept Ni,Uii and U12 constant and varied 
N2 and U22 ■ Fig. ^ plots the ratio of the total energy in the asymmetric case to the energy of the symmetric 
one as a function of these two variables. Almost all solutions are symmetry-breaking ones and a trough is 
formed near U22 — Un. The trough is an optimal region for finding symmetry-breaking solutions. 



5 Finite Interaction Range 

We have so far considered the mean-field description of a BEC mixture assuming a zero-range delta-function 
pseudopotential. It is of interest to consider the effects of a finite-range interatomic interaction on the 
ground-state structure in the two-component condensate. We introduce a pseudopotential in the form of a 
normalized Gaussian with a finite width (i.e., range) which recovers the zero-range limit result as the range 
vanishes. We search for changes in the structure of the ground state of the two-component system as the 
range of the intercomponent interactions only are varied, keeping the delta-function pseudopotential for the 
self-interactions. The results displayed below were obtained by means of direct numerical simulations, not 
the FT approximation. In ID, the intercomponent-interaction terms in Eqs. ([T]) Ui2\'4'i{x,t)\'^ are replaced 
by a nonlocal expression, 

where d is the interaction range. The Gaussian form was chosen to model a finite-range potential solely 
for its simplicity (see also Ref. for the role of a finite interaction range in attractive single-component 
BEC within this model) . This family of finite-range potentials has a constant scattering length within the 
first Born approximation. All the cases that we investigate below correspond to configurations with separated 
wavefunctions in the usual TF limit (see Sec. |3.2D . 

5.1 Weak separation {U11U22 < (Uuf) 

Here we discuss case 2 of Table 1. Fig. ^ shows the overlap region between the wavefunctions of the 
two components; this overlap region grows as the interaction range increases. Starting from relatively 
well separated phases, we end up with a complete overlap of the two components (i.e., the component 
located initially outside the narrow- width component finally penetrates it). The intercomponent interaction 
parameters are chosen as U12 = 1.02 C/n — 1.02 [/22, which places this case not far from the boundary 
between the cases of separated and overlapping phases. In other words, increasing the interaction range 



6 



forces a transition from separated phases to penetrating ones. Within the zero-range model, this would 
correspond to a transition from U12 > VU11U22 to U12 < \/UiiV22, i-e., to attenuation of the interaction. 
A simple argument justifies this conclusion. Suppose the interaction range is so large that the long-range 
potential varies only slightly over the extent of the interface between the nearly separated components 
(where the intercomponent interaction is important). Then, the interaction terms of Eq. ( |2^ ) would be 
approximately constant, producing only a shift in the total energy of the mixture, but not affecting the 
shape of the two wavefunctions. When d becomes comparable to the penetration depth (see also Ref. [^), 
phase separation is reduced. This is clearly illustrated in Fig. |^. 

Although the parameters in the ID solution are not directly relevant to an experimental situation, we 
decreased the range of the intercomponent interactions so that they are in a realistic range of values. From 
the aforementioned arguments, we see that in order to observe a difference from the zero-range case, the 
interaction range should be comparable to the penetration depth. In order for the interaction range to be 
small and yet correspond to the boundary between overlapping and segregated phases, a large disparity 
in the number of atoms in the two components is required. This situation is depicted in Fig. ^, where 
U12 — 1.01 J7ii — 1.01 C/22, but Ni N2 ■ We observe a qualitative change in the ground-state solution: 
the component that initially was at the center of the trap moves to its outskirts as the interaction range 
grows. Although c? = 0.1 y^U/mu! (bottom frame in Fig. would usually correspond to several tens of 
nanometers, we argue that one can optimize the sensitivity to d by considering the regime of parameters 
near the boundary between the penetrating and segregated phases. The possibility of manipulating the 
strength of atomic collisions is not excluded, as Feshbach resonances have been observed in BEC samples 
PH , and several other proposals in this respect have been put forward . Using such techniques one 

could vary the intercomponent scattering length in order to scan the region near detU — 0. Then, by 
comparing the measured structure of the mixture to the predictions of the two theoretical models (i.e., the 
ones with zero and finite interaction range) one could determine the effective range of interactions (which is 
a parameter in the latter model). Thus, one can probe the microscopic parameter d via a magnified effect 
such as the qualitative change of the condensate structure from to separated phase to penetrating phase. 

5.2 Strong separation {U11U22 ^ (^^12)^) 

In the preceding section we demonstrated the effective attenuation of the mean-field repulsive interaction 
between two components of a BEC due to an increase in the range of the interaction for the case when 
UiiU22 < (£^12)^- Now we turn to the case when the parameters of the BEC mixture are far from threshold 
for the onset of the penetrating phase. In Fig. ^ is for parameters U12 = TUn = 7U22, hence the interface 
between the two components is very sharp. Since now it is much easier to match the range of interactions 
with the penetration depth, one might expect that effects similar to those described above will appear at 
even smaller values of d. However, this is not the case. The mutual repulsion of the components remains very 
strong even if reduced by a finite interaction range. As the interaction range increases, the two components 
tend to move apart, yielding two completely separated phases. 



6 Summary and Conclusion 

We have presented a detailed classification of stable solutions for binary mixtures of dilute atomic conden- 
sates. The analysis is particularly simple within the Thomas-Fermi approximation. Within this approx- 
imation one can distinguish two general classes of ground state for two component condensate mixtures: 
unseparated ones with an overlap region (both component wavefunctions are simultaneously nonzero) and 
separated ones not containing an overlap region (except for the tail penetration). The latter contains 
also solutions that do not posses the symmetry of the trapping potential. Components are separated if 
detU = U11U22 < U12 < and they overlap if detU > 0. The predictions from the TF approximation 
become ambiguous in the region of parameters where phase separated solutions which break the symmetry 
of the trap are energy-degenerate with the phase-separated solutions preserving the symmetry. In this case, 
it is crucial to include the contribution of the kinetic energy operator and of the mutual interaction energy 
in determining the structure of the ground state geometry, which tend to favour the asymmetric solutions. 
The physical reason for this is probably a smaller interface region in the asymmetric solutions. It is these 
interface regions which contribute most to both the kinetic and mutual interaction energies. We have carried 
out our numerical calculations in ID, but our conclusions should be valid in two and three dimensions as well. 
The condensates, if not overlapping, should have a propensity towards states with minimal interface surface 
area. In the det U < case, with kinetic energy included, we further recognize a weak separation regime 
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(det U < 0) and strong separation regime (det [7^0) when a penetration depth goes to zero. Contribution 
of the kinetic energy to the total energy increases with decreasing det U (for negative det U) , in spite of the 
decreasing interface size (see Fig. (^). This is due to the increasing importance of the U12 term that gives 
the cross-interaction energy of atoms from different components. 

Since the size of the overlap region can be very small (smaller than a single condensate healing length), 
it is of interest to investigate the possible impact of a small, but non-zero, interaction range in a binary 
condensate. We have developed a model of finite range potential by introducing a pseudopotential in the 
form of a normalized Gaussian with finite width. We identified two distinctly different cases. In one case, that 
of significant overlap, the finite range tends to increase the penetration over the delta function interaction, 
and in the other, that of strong separation, the finite range leads to a trough between the two condensates. 
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by the US-Israel Binational Science Foundation and the James Franck Binational German-Israel Program 
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Figure 1: Two types of solutions with partial overlap. Plotted are densities for the first and the second 
component as solid and dashed curves. Panels (a) and (b) depict have Un : U12 : U22 ~ 1-2 : 0.9 : 0.8 and 
Uii '■ U12 : U22 = 1-1 : 0.9 : 1 respectively. In both panels A^i = A^2- The oscillator length unit is given by 
y^h/muj. Density distributions in all figures are normalized to unity. 
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Figure 2: Two types of solutions without overlap: symmetric and asymmetric cases. Plotted are densities 
for the first and second components as solid and dashed curves. For both panels, the ratio of the scattering 
lengths is Uu : U12 : U22 = 1 : 1-52 : 1.01 and Ni = N2. The ratio of the energy of the asymmetric case to 
the symmetric one is 0.8. 
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Figure 3: Total energy vs. the dimensionless quantity AeXlJ j {UwV 22) for the symmetric and asymmetric 
cases (numerical simulation) and the TF predicion. Here N\ = N2, Un and U22 are kept constant (the same 
as in Fig. whereas U12 is varied. For arguments smaller than the three curves are indistinguishable. 




Figure 4: Ratio of the total energy of the asymmetric case to the corresponding value in the symmetric 
phase vs. U22 and N2 scaled by certain initial values. The initial values are the same as in Fig. |^ (i.e., 
Uii : U12 : C/22 = 1 : 1-52 : 1.01 and Ni = N2). Note the deep valley near Un = U22. 
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Figure 5: Dependence of the density distributions of a two-component BEC in weakly segregated phases on 
the range of the intercomponent interaction d. Here U12 = 1.02 i/n = 1.02 C/22 and N2/N1 = 0.15. The 
effect of attenuation of the interaction with growth of the range is illustrated for three values of d. 
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Figure 6: Interaction-range induced change of the ground-state density distribution of a two-component 
BEC in weakly segregated phases, but with a large disparity in the number of atoms in the components. 
The parameters are U12 = 1.01 Un = 1.01 U22 and N2IN1 = 0.02. 
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Figure 7: Enhanced mutual repulsion as an effect of a finite interaction range between strongly separated 
phases. Here Uu : U12 : U22 = 1 ■ 7 : 1 and N2/N1 = 8. 
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